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We develop the concept of scattering matrix and we use it to perform stable numerical calcu- 
lations of resonant tunneling of electrons through a multiple potential barrier in a semiconductor 
heterostructure. Electrons move in two external nonperturbative electric fields: constant and os- 
cillating in time. We apply our algorithm for different strengths and spatial configurations of the 
fields. 
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I. INTRODUCTION 



The aim of this paper is to present a numerically stable 
algorithm for investigations of nonrelativistic quantum 
processes occurring in arbitrary space-dependent scalar 
potential and a time- and space-dependent vector poten- 
tial. Vector potential is periodic in time and describes 
a laser field. Such conditions are met for example in 
semiconductor nanostructures [lHal (like quantum wires 
or wells), carbon nanotubes |9lTl0j or in surface physics 
[Tll - [T5j . To make our presentation as clear as possible we 
shall restrict ourself to the one-space-dimensional case, 
although extension of this algorithm to two and three- 
space dimensional systems, also with magnetic field ac- 
counted for, is possible (see, e.g. [lf|) . We shall ap- 
ply our method to investigation of the tunneling process 
and its dependence on relative phases of multi-chromatic 
laser pulses (multi-color processes have been considered 
for instance in [TBI. Kt\). 

Multiple barrier, field-assisted resonant tunneling is an 
interesting problem because it provides an insight into 
the physics of nanostructure quantum systems and be- 
cause it is a fundamental effect to use in a wide vari- 
ety of technological applications. As to the latter, it is 
enough to mention all sorts of detectors and generators of 
microwave radiation based on double barrier structures 
with external electric field added; for more examples, 
see [lS - 23] . Here we analyze resonant tunneling through 
semiconductor structures in the presence of both oscil- 
lating and constant in time external electric fields. The 
fields are assumed to be nonperturbative. We assume 
that single-particle states of electrons in heterostructures 
are well approximated by the so-called envelope func- 
tion [2^, [23| . Effects of sharp interfaces between different 
semiconductors are accounted for by boundary conditions 
satisfied by the envelope wavefunction, i.e., by the conti- 
nuity of both the envelope wavefunction and the proba- 
bility current at the interfaces (see, for example, [24143 1| ). 
Scalar potential V(x) is assumed to be constant in time 
but it can be of any shape. The same conditions hold 
for space-dependent effective mass m(x). Vector poten- 
tial A(x, t) describes a laser field and thus it is space- 
dependent and oscillates in time. In our approach to 



numerical computations, one-dimensional space is sliced 
into small intervals where m(x), V(x) and A(x,t) are 
space-independent. For arbitrary space-dependent func- 
tions m(x), V(x) and A{x,t) such a procedure is justi- 
fied provided that the widths of these intervals are suf- 
ficiently small. We develop below a general numerical 
scheme which permits to evaluate transition and reflec- 
tion probabilities for electrons moving in the system de- 
scribed above. 

This paper is organized as follows. In Sec. II the 
most general solution of the Schrodinger equation is in- 
troduced. The transfer-matrix method and matching 
conditions are analyzed in Sec. Ill, whereas reflection 
and transition probabilities are introduced in Sec. IV. 
These probabilities must sum up to 1, which puts a very 
strong check for the accuracy of numerical calculations. 
The most important part of this paper, i.e. the concept 
of the scattering-matrix method, is discussed in the next 
section, where it is shown why the scattering-matrix al- 
gorithm has to be introduced, instead of a much simpler 
transfer-matrix algorithm. Numerical illustrations of the 
applicability of this algorithm are presented in Sections 
VI and VII, and are followed by short conclusions. 



II. SOLUTION OF THE SCHRODINGER 
EQUATION 

Let us start with one-dimensional Schrodinger equa- 
tion of the form I20H2H, 



id t yb{x, t) = \\ (-d x - eA(x, tj) — ^- (-d x - eA{x, tj) 
L2 V 1 / m(x) V 1 / 



+V{x) ip(x,t). 



(1) 



Space-dependent mass m(x), scalar potential V(x) and 
vector potential A(x, t) are spatially constant in finite 
intervals. Their values in any interval (xj_i,Xi) will be 
denoted as mi, Vi and Ai(t). An example of such a struc- 
ture is presented in Fig. [1] We require also that the 
function A(x,t) is periodic in time, that is 



A{x,t + T) = A(x,t), 



(2) 
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where T = 2-k/uj and uj is the frequency of the oscillating 
in time electric field. Defining in a standard way the 
probability density p(x,t), 

p(x,t) = \t/j(x,t)\ 2 , (3) 
and the probability current j(x,t), 

j(x,t) = \r{x,t)-^(\d x -eA{x,t))^{x,t) (4) 
2 m[x) V i / 

f 



\${x,t)—— 
2 m(x) 



(t* 



eA(x,t))ip(x,t) 



we show using Eq. ([TJ that the conservation of proba- 
bility condition is satisfied. Indeed, assuming the above 
definitions, we get the continuity equation, 



d t p{x,t) + d x j(x,t) = 0. 



(5) 



Space dependence of mass in Eq. (JTJ forces one to im- 
pose non-standard continuity conditions on any solution 
of this equation. It is now the wavefunction ij)(x,t) and 
the quantity 



1 (ld x - eA{x,tj)i>{x,t) 



m(x) 



(6) 



that have to be continuous at points of discontinuity of 
mass m(x) and both potentials V(x) and A(x,i) (25142^ . 
Before passing to a general solution ij){x,t) of Eq. flTJ 
in any given interval (x^-i, Xi), which we shall denote as 
tpi(x,t), let us note that due to time periodicity of the 
Hamiltonian, ipi(x, t) can be chosen such that the Floquet 
condition, 



(7) 



is satisfied, where E is the so-called quasienergy. A gen- 
eral solution tpi(x, t) of Eq. ( []]) in any interval Xi) 
takes then the following form [29|, [30] , 



ipi(x,t) = exp(-i(E + Muj)t)Y^ (8) 

M= — OO <T=± 

oo 

^ C° n B M -n (crPiN) exp (iap iN x) , 



N=—oo 



where Cf N are arbitrary complex numbers to be deter- 
mined and 



PiN = y/2mi(E + Nu-Vi-Ui), 



(9) 



with Ui — e 2 (Aj (t)) / 2nii being the ponderomotive en- 
ergy, where (A 2 (t)) means the time-average of A\ (t) over 
the laser-field oscillation. Components for which are 
purely imaginary are called closed channels. These chan- 
nels are not observed for a particle in initial or final 
states, but they have to be taken into account in order to 
satisfy the unitary condition of the time evolution. In a 
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FIG. 1: Generic shapes of space-dependent superlattice po- 
tential V(x), effective mass m(x), and oscillating in time laser 
field A(x,t). 



general case, the Bm-n{<^P%n) functions are components 
of the following Fourier expansion, 

oo 

exp(i$? N (t)) = ex P {-iMwt)B M ^ N (<jp m ) (10) 

M=-oo 

provided that the vector potential A(x, t) is periodic in 
time. Functions ®f N {t) are defined as follows: 

n N {t) = / —Mt) PlN - —(A 2 (t) - (A 2 ( t ))) d*. 

It is easily seen from the above equation that the 
Bm-n^Pin) functions depend on the form of the vector 
potential A(x,t), that is on the laser field applied. 



III. MATCHING CONDITIONS AND 
TRANSFER MATRIX 

Continuity conditions discussed above and applied to a 
general solution (|SJ of the Schrodinger equation ([J} lead 
to an infinite chain of equations connecting constants Cf N 
in the neighboring domains. These matching conditions 
can be written in the matrix form, 

B(i - l.ii.iJCU = B(i, Xj_i)Cj, (12) 

where Cf N = [C i ]n are the components of the columns 
Cf. The matrices B(i,x) and Ci are defined as follows, 

B(i x) ( K {hx) B ~ {hx) \ C -( C ?\ (13) 

The elements of B(i, x) can be computed in the following 
way. 
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For an arbitrary function A(x, t), periodic in time with 
the period T, 



A(x,t) =A(x,t + T) 



(14) 



we have 



A(x,t) = b n (x) exp (— incot), (15) 



where uj — 2tt/T. In the interval (xi-\,Xi) coefficients 
b n {x) assume constant values, which we shall denote as 
&j >n . Using the condition of the continuity of the wave- 
function ipi(x,t) at the point a;,_i, we compute the ele- 
ments of the matrices B + and B~ , 

B ± (i,x) M ,N = B M -N(±Pi,N)exp(±ipi,Nx)- (16) 

On the other hand elements of the B' matrix can be 
evaluated by substituting a general solution ((5|) to the 
expression ([6]) and applying the continuity condition to 
it at Xi-\. After some algebraic manipulations we obtain 
the following equation, 

f V exp (-UE + Mu))t) 

1 1 M=-oo 
oo 

x E E C°_ lN B M -N(vpi-i,N)aPi-i,N 

ty=±l N=~oa 

oo 

x exp (\api-i.NXi-i) — exp (— i(E + Mui)t) 

M=-oo 

oo 

-iV-n(oPi-l,iv) 

a=±l N,n= — oo 

x expiiapi-i^Xi-xfj 

1 

= — ( V exp (-i(E + Mute) 

M=—oo 
oo 

X E E C i,NBM-N(<JPt,N)vPi,N 
a=±l N=-oo 

oo 

exp(iapi t pfXi-x) — ^ exp (-i(E + Mu)t) 

M= — oo 

oo 

x E E eb^ n Cl N B M 



er=±l N,n— — oo 

x exp {\api, N Xi-i) 



(17) 



Suspending the summation over M on both sides of the 
above equation, we finally get the expression for the B'- 
matrices, 

B' ± (i 1 x) M ,N = ±—B M -N(Pi,N)Pi,Nexp(±ipi. N x) 
mi 



- — V, e.b ln B M -N-n{±Pi,N) 



n— — oo 



In this way we obtain a set of equations for vectors C^, 
C i = B i C i - 1 , (19) 

where 

B t = [B&Xi-i^Bii - 1, £,_!). (20) 

These relations allow to connect a solution in a given 
domain Xi-\ < x < Xi with an analogous solution in any 
other domain Xj-\ < x < Xj, 



— BjBj-i, . . . , Bi + \Ci — TjiCi 



(21) 



where Tji is the so-called transfer matrix [2JJ, |27J, |29, 31. 



IV. REFLECTION AND TRANSITION 
PROBABILITIES 

It is clear now that on the basis of Eq. (|2"T1) we can con- 
nect solutions in the boundary domains (— oo,ieo) and 
(xl_i,oo). Values of mass m(x), scalar potential V(x) 
and vector potential A{x, t) in these domains will be de- 
noted as too, Vb, Ao(t) and tol, Vl, Ai,(t), respectively. 
We can then write down solutions of ([T]) for each of these 
domains. These solutions represent incident (V'mc), re- 
flected (?/> re f) and transmitted (i/'tr) waves, and take the 
following form, 

oo 

vpi nc (x,t) = exp (— iEt) exp (— iMuit) 



M=-c 



x B M (po)exp(ip x), 



(22) 



t[> TC f(x,t) = C n exp (— iEt) exp (— iMuit) 

N,M=-oo 

x B M -N(~PN)exp(-ip N x), (23) 

oo 

ip tr (x,t) = Yj C£ N exp (— iEt) exp (— iMuit) 

M=-oo 

x B M -N{qN)exp(iq N x), (24) 

where 



PN = ^2m (E + Nuj-Vo-U ), 

q N = ^2m L {E + NLo~V L ~U L ). (25) 

Constants Cq N and C\ N will be denoted from now on 
as Rn and Tjy, respectively. Using continuity conditions 
for functions defined above, we get the probability conser- 
vation equation for reflection and transition amplitudes, 
Rjv and T^r, 



x exp (±ip iiN x) 



(18) 



E -iR-i 2 _ 



E 



mo9JV| T^| 2 = l, (26) 
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where summations are over such N for which p^ and 
§jv are real, i.e., over the open channels. This equation 
permits us to interpret 



THE SCATTERING MATRIX 



Pr(A0 = — \Rn\ 2 
Po 



and 



Pt(N) 



1 ijv 



(27) 



(28) 



m L p 

as reflection and transition probabilities for a tunneling 
process in which absorption (N > 0) or emission (N < 0) 
of energy Noj by electrons occurred [H, [29[ . In case of a 
monochromatic laser field this process can be interpreted 
as absorption or emission of N photons from the laser 
field. 

The unitary condition (f2l)f can be also interpreted as 
the conservation of electric charge. To this end, let us de- 
fine the quantities proportional to the density of electric 
currents, 

Po 



J'mc — 
J rcf — 

Jtr = 



too 

E -iR-i 2 ' 

too 



(29) 
(30) 



E — i T -i 2 - ( 31 ) 

N^Ntr 

Then Eq. ([26)) adopts the form of the first Kirchhoff low, 

Jinc = Jiei + Jtr- (32) 

Using (f2"Tj) we can calculate constants C^ N = Rat and 
since 



Tn appearing in equations (|22|) - (|24|) . Indeed, 



C L = TCo, (33) 

where transfer matrix T = Tlo, and because T, Co and 
Cl adopt the following block forms, 



T 



T++ 

r-+ 



T 



— ,0> 



we arrive at 



T 




r- + c' 



R 



,C L 

R, 
R, 



Thus, after some algebraic manipulations, 



where R and T denote the columns of Rjv i T 

[Co"]jv = So,N' 

we have, 

R = 
T = 



N 



,(34) 



(35) 
and 



-(T- 
(T+ 4 



)- x r 

T + -(T 



°o ■ 



o ■ 



(36) 



which allows us to determine the quantities Rjv and 
Tjv for a given transfer matrix T . For open channels, 
these quantities are the amplitudes of reflection (Rat) and 
transition (Tjv) probabilities, from which one can com- 
pute reflection and transition probabilities using equa- 
tions (|2~T|) and (|2"8j) . In all our numerical illustrations, 
condition (|26p is satisfied with an error smaller than 



Wc note from equations (fT6|) and (fT8|) that each of the 
Bi matrices that constitute the transfer matrix Tji con- 
tain elements exp(±ip,* ) jva;i) that depend on the Xi coor- 
dinates at which the discontinuities appear. For closed 
channels, that is when the momenta are purely imag- 
inary, these numbers are real and may assume arbitrary 
values, very large or very small, depending again on the 
Xi coordinates. Number of the Bi matrices is equal to 
the number of discontinuity points, that is it depends on 
how we divide the space into short intervals in order to 
make our potential tractable by our algorithm. It may 
therefore turn out that in order to compute the transfer 
matrix Tji , we have to multiply a large number of the Bi 
matrices, each containing both very small and very large 
numbers. It is clear that such a procedure is numerically 
unstable. We have to find a way to modify our method 
of calculations in order to compute the elements of each 
Bi matrix at the same point x = independently of 
where the 'real' Xi is. This would eliminate "dangerous" 
exTp(±ipi^]\[Xi) elements (turning them to 1), however at 
the cost of appearing somewhere else. We shall see later 
that these 'left-overs' of the shift into x = appear only 
as differences Xj+i — Xi and therefore do not cause any 
harmful side-effects. We shall see now that such a mod- 
ification is possible and the price we pay for it is worth 
the effort. 

It follows from Eq. (l2~Tj) that in the neighboring do- 
mains, (£&j_2, Xi—i) and (xi-i,xi), we have, 



Cj — TiA~\Ci-\. 



(37) 



Although the elements of the transfer matrix 7i,i-i have 
been computed from the continuity conditions at point 
ajj-i, one can compute them at any other point, for ex- 
ample x — 0. To this end, let us notice what follows 
from the form of the solution (|8|). Translation of the sys- 
tem by a certain distance 8 along the x-axis causes only 
multiplication of each member of the sum over N in ((5J) 
by a constant exp (iapnyS). These constants can be in- 
cluded in coefficients Cf N . In this way we get a new set 
of constants which we shall denote as C% 



N- 



Cf N = exp (iapiN^C^ 



A'- 



(38) 



We shall interpret these constants as coefficients in so- 
lution ([8]), given by the continuity conditions at point 
Xi—i — 5. Eq. (|38p written in the matrix form becomes, 



where 



and 



-14 



Q = Pi(6)Ci 



Vi(s) = (pm o 

V Pr(5) J ' 



(39) 
(40) 

(41) 



5 



o 

3.2 



X, 



V 2 {x 2 -x 1 



'2,1 

-o- 

3> -I 



'1,0 



PlfXi-Xo) 



FIG. 2: Diagrammatic representation of Eq. (|50[) . Circles 
represent points of discontinuity {xj} and matrices {TJ+i^ }, 
whereas lines represent 'free-propagators' {Vj+i(xj+i — Sj)}. 
It is important to notice that all matrices {Tj+i A are calcu- 
lated at x = 0, which prevents the development of numerical 
overflows. 
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In the equation above Pf{6) is a diagonal matrix, 

[P?($)]nn> = <W exp (iap iN 5), (42) 
whereas Cr and Cr are the columns of the constants Cfi 



Cf N and [Cf] 



>N 



N 



and Cf N respectively, that is [C^]n 
Cf N . It follows from the form of the matrix Vi (6) that 
the following relations are satisfied: 



Pi(Si)-Pi(S 2 )=V i (6 1 + 6 2 ) 



(43) 



(44) 



Let us notice also that translation of the system defined 
above modifies the transfer matrix Tm-i. We have 



Pi 1 Ci — Ci — Ti t i- 

= Ti.i-iPi^Vi-ii^d-!. 



-lCj-i 
,-1 



thus 



= P t (6)%,i-iP^\(6)C t - 1 , 
and we can write it down as 

Ci = 77,j_lVi_i, 

where 



(45) 



(46) 



(47) 



(48) 



Matrix elements denoted with the tilde symbol refer to 
the translated system. Using the method defined above 
and the relation (|2~T1) . we can connect now the solu- 
tion in the domain (— oo,xo) with the solution in any 
other domain {xi-\, Xi). In this way the elements of the 
transfer matrix, which have been computed until now at 
the points of discontinuity xq . . . X{-\, are computed now 
each time at the same point x = 0. Let us illustrate this 
method for a special case of i = 3 

C 3 = % a Ti,xTxfiC x =Pz 1 {x 2 )^ 2 P 2 {x 2 )P 2 - 1 {x 1 ) 
x T 2 " 1 Px(xx)P 1 - 1 (x Q )T 1 Po(x )C 

= P3 1 (^2)^ 2 P2(X2-X 1 ) 

x T^Pxixi - x )7?. Q Po(xo)C . (49) 



FIG. 3: Schematic representation of the idea of the transfer 
matrix and the scattering matrix. For the transfer matrix the 
incoming channels are Cf and C~ , and the outgoing channels 
are and CJ . For the scattering matrix Cf and C~ are 
the incoming channels with the remaining two considered as 
the outgoing ones. 



Equation (|4"9l connects constants Cq and C3 using the 
matrices Tj j-i all computed at x — independently 
of j, and diagonal matrices Vj(5j), given by the rela- 



tions (j40j and (|42|). where Sj = (xj 



Edge ma- 



trices Vq(xq) and P^ 1 (x 2 ) in the equation (|4"9")l can be 
omitted while computing the transmission and reflection 
probability amplitudes since their only role is to multi- 
ply the amplitudes by phase quotients which disappear 
while computing the probabilities. Although these matri- 
ces lead to significant modifications of the closed channels 
in the domains of x < xq and x > X3 in this particular 
case, these channels do not influence the reflection and 
transition amplitudes. Transmission and reflection prob- 
abilities can thus be computed using a modified transfer 
matrix, 



%° = T*,P 2 (x 2 - x x )TtxVx{xx ~ *o)7? . 



(50) 



The matrices 7~°i_i are equal to the matrices Pi in Eq. 
(|2U)l calculated however for Xi-x = 0. This fact speeds 
up numerical calculations since now matrix P(i,x = 0) 
in Eq. (|20"|) have to be inverted only once. Further on we 
shall omit the superscript in T and the tilde over C in 
order to simplify notation. Diagrammatic representation 
of the equation above is shown in Fig. [5J 

The method presented above is still numerically unsta- 
ble. The reason for this instability lies in the existence of 
large numerical values of elements of the P^(5) matrix 
for imaginary momenta PiN- In other words, for 



Ci = 



V, 



= 71 



T~ + T 

'1,1 — 1 ' a 



(-1 



xCi-X 

c~_ 



(51) 



the source of numerical instabilities are matrix elements 
7~i~Zi that contain large numbers. There is however a 
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chance for improving the stability, if only its reverse will 

be used, (7~~Zi) ■ It appears that it is possible pro- 
vided that in our numerical algorithm only the so-called 
scattering matrix will be applied. For this reason we will 
show below how to compute the scattering matrix, Sj t i, 
using only elements of the transfer matrix, Tj,i- For the 
transfer matrix T 3 I j we have, 



Tj : iCi — Cj 
Thus, 



c: 



T++ T + - 

3,1 , 3,1 

Tr.+ r~- 



c; 



.(52) 



C7+ = 7-++C+ + 7-+-C7 



C; 



7j,i c i 



(53) 



On the basis of (|53|) we now want to compute the el- 
ements of the Sj i matrix. This matrix is supposed to 

—i -t- 

connect the coefficients C, and Cj~ in the following way 
(for the graphical illustration, see Fig. [3]), 



CJ 



%. — 



C t ) l S j/ S 3,i ) ( C 3 



(54) 



Using the set of linear equations (|53|) , we easily compute 
the coefficients C~ and on the left-hand side of equa- 
tion (|54|) , as functions of the coefficients CJ and Cf . Wc 
get then the following relations, 

c i = WrryKCj -Tr+C+), 

C j = (Tj.i ) 1 T ji + )C+ 

+ T+r( T rr)-^Cr. (55) 



Finally we compute the elements of the matrix S 



3,1 ■ 



s t + = -VZi-)- l r: 



3,1- 



(7? 



3,1 3,i v 



3,i v 3,1- ) 



(56) 



As expected, the matrix Sj.i contains only numerically 
stable elements 

It follows from Eq. (f2Tj) that the transfer matrix Tj.i 
can be written as the product of two transfer matrices, 
Tj,k and Tk.i (i < k < j), 

Tj.i — Tj,kTk,i, (57) 

where matrices Tj, k and 7~k,i are defined as follows, 

Ck — Tk,iCi, 

C 3 = r j>k C k . (58) 

Applying the method presented above, for each of the 
transfer matrices 7},fc and Tk,i we can construct a scat- 
tering matrix, S ^k and Sk i respectively. Elements of 
the scattering matrix Sj i can be computed using only 
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FIG. 4: Tunneling process considered in this paper. Param- 
eters for the triple barrier are: Vo = 237meV, and the effec- 
tive masses mGaAs = 0.0667m e and mGa x Ai 1 „ x As = 0.0918m c , 
where m e is the electron rest mass. The widths of the barriers 
b and wells a can change. 




-100 



-50 50 

x (in angstrems) 



FIG. 5: (Color online) Plot of the potential V(x) + Fx for 
a = 40A, b = 20A, and F = -0.23 x l(T 4 (in atomic units). 
Other parameters are the same as in Fig. [4] 



elements of the matrices Sj t k and Sk t i- Using the nota- 
tion above, we obtain the following expressions for the 
elements of the Sjj matrix, 



C++ _c+H 
—°k,i 



S k ,i (1 - Sj~, k Sk,i ) S?,kSk,ti 

S j,i =,S j,k i 1 ~ S j,k S kA ) S k,t 



~ S k,i i 1 ~ S £k S k,i ) S £k > 



Sj,k + Sj,k + Sk,i C 1 _ Sj,k $k,i ) $tk 



(59) 



It is clear from the above that the Sji matrix is not 
merely a product of two matrices Sj t k and Sk,i, but rather 
a complicated nonlinear composition of them. It is im- 
portant however to note that despite its evident complex- 
ity, such a construction of the scattering matrix is numer- 
ically stable, as opposed to the transfer matrix method 
which fails if a system with a large number of discon- 
tinuity points Xi is considered. Stability of such an al- 
gorithm has been proven in our numerical investigations 
by checking that the condition (|26|) is satisfied with an 
error smaller than 1CU 14 . Such an accuracy can never be 
achieved for systems with a large number of discontinuity 
points if the transfer matrix is applied. 
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100 200 300 

E (meV) 

FIG. 6: (Color online) Transmission probabilities for the po- 
tential shown in Fig. [S] Numbers in the insert frame (repre- 
senting enlarged part of the main frame) indicate the number 
of equally spaced discontinuity points introduced in our nu- 
merical algorithm. We see that 15 points do not give correct 
results and that the convergence is reached with more than 
100 points. 



VI. RESONANT TUNNELING 

We shall consider now the tunneling phenomenon 
through a semiconductor heterostructure presented in 
Fig. |U In the beginning, let us assume that electrons 
interact only with a constant electric field. Hence, the 
time-independent potential is of the form V(x) + Fx, 
in which V(x) represents the semiconductor heterostruc- 
ture potential (Fig. |4} and F is the electric-field strength. 
The plot of this potential is presented in Fig. [SJ where 
a = 40A, b = 20A, and F = -0.23 x 10" 4 (in atomic 
units). Applying the theory developed above, we calcu- 
lated reflection and transmission probabilities (see Fig. 
[5]) discretizing the potential with 15, 141, and 281 points, 
as indicated in one of the frames. We observe that, in 
order to get convergence, one has to introduce at least 
one hundred discontinuity points. There is no noticable 
difference between the results obtained for 141 and 281 
such points. 

Next, let us analyze transmission of electrons through 
the triple barrier of Fig. g] with a = 70A and b = 20A 
and with the 221 discretization points. Now we apply a 
constant electric field and the monochromatic laser field 
of frequency w = 70meV and intensity such that its pon- 
deromotive energy divided by the laser photon energy 
equals 10~ 4 . Without external fields, the resonant ener- 
gies are grouped in doublets in which the lower-energy 
resonance corresponds to the antisymmetric resonance 
state, and the upper-energy resonance to the symmetric 
one. With the laser field switched on this structure does 
not change very much provided that the frequency is off- 
resonance with respect to the already existing resonance 
states of the triple barrier, and the intensity is not too 
large, as it is presented in Figs. [7] and [5] The pattern 
changes significantly if a constant electric field is applied. 




E (meV) 

FIG. 7: (Color online) Transmission probabilities for semicon- 
ductor triple barrier with a = 70A, b = 20A. Intensity of the 
laser field is such that the ratio of ponderomotive energy to 
photon energy is U p /uj = 10~ 4 with to = 70meV and we have 
three electric- field strengths (in atomic units), as indicated in 
the figure. As expected, transition probability distributions 
for ±F [blue (dash-dash) and red (dash-dot) lines] are shifted 
in energy by |f |(36 + 2a). Computations were performed for 
221 discretization points. 



We observe that with an increasing strength of the elec- 
tric field the low-energy transmission resonances from a 
given doublet gradually disappear and we are left with 
a single transmission resonance, which for even stronger 
electric fields also dies out. This means that by proper 
tunning the strength of a constant electric field one can 
selectively transmit electrons of some particular energies. 



x10" 5 
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FIG. 8: (Color online) Color map of the total transmission 
probability in the plane of the incident electron energy E and 
the electric-field strength F. As expected, for the vanishing 
electric field resonances in the considered in Fig. 0triple bar- 
rier structure show up in doublets. However, for sufficiently 
strong electric field the lower-energy resonance from a doublet 
disappears. 
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FIG. 9: (Color online) Total transmission probabilities for 
the triple barrier structure presented in Fig. 0]with a = 70A 
and b — 20A, and for the bichromatic field [Eq. (|60[) ] with 
ui — 70meV. The continuous (blue) line is for tp = 0, dash- 
dash (red) line for tp = n/2, whereas dash-dot (black) line 
for tp = 7T. Frames correspond to different laser field in- 
tensities characterized by the dimensionless parameter £ = 
\e\Eo(x) / (2\/ hrricUi 3 ) (m c is the electron rest mass and Eo(x) 
is considered to be constant in the whole space): (a) £ = 0.1, 
(b) £ = 0.5, (c) ^ = 1 and (d) £ = 2. 



VII. PHASE CONTROL OF TUNNELING 

Special features of barrier problems stem from the in- 
teraction of waves reflected from or transmitted through 
potential jumps. When the interference of reflected waves 
is in phase, transmission becomes minimal. But when the 
interference of reflected waves is out of phase (i.e., it is 
destructive) the incident wave resonantly penetrates ei- 
ther by tunneling through or passing above the barrier 
structure. If the process occurs in a monochromatic laser 
field the destructive or constructive interferences between 
reflected and transmitted waves are present also for dif- 
ferent Fourier components of the electron wavefunction. 
This leads for example to opening or closing gaps in the 
band structure P-d, H3] or formation of multiple-plateau 
structures in the high-order harmonic spectrum |33l . [34| . 
It gets even more complicated if multichromatic laser 
fields or short laser pulses are applied. In the first case 
the interference discussed above can be controlled by rel- 
ative phases of harmonics present in the multichromatic 
fields, whereas in the second case the resonance trans- 
mission can be modified by the carrier-envelope phase. 

As an example let us consider a bichromatic laser field. 
Let the electric field be of the form 

E(x,t) = E {x)[sm(tut) - sm(2ut + cp)], (60) 

where Eq(x) is in general a space-dependent amplitude 
of the laser field. In Figs. H] and [TU] the laser-modified 
total transmission probabilities through a triple-barrier 



FIG. 10: (Color online) The same is in Fig. [9] but with space- 
dependent intensity of a laser field. Now, Eo(x) = Eo within 
the triple barrier structure, Eo(x) = Eo/2 within the edge 
barriers and outside. The electric field strength Eq is deter- 
mined by a dimensionless parameter £ = |e|i?o/ (2y/ hm c ui 3 ) , 
with the same numerical values as in Fig. [9] 



structure are presented for three different relative phases 
tp. Fig. IH1 corresponds to the situation in which the laser 
field acts in the whole space, whereas Fig. [TU] illustrates 
the action of the laser field concentrated within the struc- 
ture, hence incident, reflected and transmitted electrons 
are free. Apart from a significant dependence of these 
probabilities on the relative phase we observe also that 
in the second case and for higher intensities considered 
the transmission probabilities are smaller. It is because 
the electrons have to traverse an extra 'potential barrier' 
created by the ponderomotive energy of a laser field in 
order to transmit through barriers. This finding opens 
up a possibility to create tunneling barrier structures by 
laser fields modulated in space. This can be investigated 
numerically by applying the algorithm developed in this 
paper. 

As a second example of the phase control let us con- 
sider transmission of electrons through a triple barrier 
structure in the presence of both a constant electric field 
and a train of very short laser pulses. It is well-known 
from atomic and molecular physics that the ionization 
process can be significantly modified by the so-called 
carrier-envelope phase if a single pulse contains only few 
oscillations. In order to investigate this phenomenon 
for electron transmission let us assume that the train 
of pulses is build from a single pulse (defined for times 
«S t T p ) of the form 



E(x, t) = E (x)f(t) sm(ujt + tp) 
where the envelop function f(t) equals 



(61) 



f(t) = exp 
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FIG. 11: (Color online) Total transmission probabilities for 
the triple barrier structure presented in Fig. [4] with a = 70A, 
b — 20A and Vo = 237meV. The train of laser pulses [Eqs. 
([61) and ((62j)] with uj = 70meV, T p = 26tt/cj and a p = T p /140 
(one-cycle pulse) is defined by the space-dependent electric 
field such that Eo(x) = Eq within the triple barrier struc- 
ture, Eq(x) = -Bo/2 within the edge barriers and outside. 
The laser field intensity is characterized by the dimension- 
less parameter f = \e\Ea/{2\/hm c Lo' A ) (m c is the electron rest 
mass), whereas the constant electric field strength F is de- 
termined by the parameter r\ = |e|-F(36 + 2a) /Vo- In this 
illustration £ = 0.2 and r\ = 0.1. The continuous (blue) line 
is for (p = 7r/2, dash-dash (red) line for <p = 7r/4, whereas 
dash-dot (black) line for tp = 0. 



and the constant in time is chosen such that 

/ " E(x,t)dt = 0. (63) 
Jo 

The carrier-envelope phase tp can change from to 2ir. 

In Fig. [Tl]we present transmission probabilities for the 
electrons impinging from the right on the triple barrier 
structure shown in Fig. [5] (however, with different val- 
ues for a and F). Without the action of the laser pulse 
the transmission is forbidden for energies smaller than 
approximately lOOmeV, whereas for larger energies elec- 
trons can tunnel resonantly. The presence of the laser 



field modifies these conditions and they get similar to 
those met in the photoemission from solid surfaces or 
the ionization of atoms or molecules. For the latter it is 
well-known that the carrier-envelope phase substantially 
modifies ionization probabilities [35l [36| . The results 
presented here also confirm this effect for the tunneling 
phenomena. Transmission probabilities (hence, also pho- 
tocurrents emitted from the surface) can then be changed 
by the carrier-envelope phase even by two orders of mag- 
nitude. 

VIII. CONCLUSIONS 

As mentioned above, our algorithm is convergent pro- 
vided that a sufficient number of discretization points is 
introduced. For systems considered here, this number 
should not be smaller than 100. If the laser field is very 
weak, this does not create significant numerical problems, 
except that calculations become longer. However, when 
the laser field is sufficiently intense, the algorithm based 
on the transfer matrix is unstable. This instability is due 
to the existence of closed channels, which introduce into 
numerical calculations very small and very large numbers 
at the same time. Augmenting precisions significantly 
slows down the calculation and does not diminish the 
problem. We have found that it is possible to make this 
algorithm numerically stable by just applying nonlinear 
matrix transformations, without introducing higher pre- 
cisions. 

Illustrations presented in this paper show that tun- 
neling of electrons through multi-barrier semiconductor 
structures can be changed significantly by applied non- 
perturbative electric fields: oscillating in time or con- 
stant. The efficiency of the algorithm presented in this 
contribution opens up the possibility of investigating sur- 
face phenomena (like photoemission or high-order har- 
monic generation) in the presence of more realistic laser 
pulses that gradually decrease within solids and extend 
on a mesoscopic scale in vacuum. These problems are 
under investigations. 
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